**# 0. settings, user and directories
clear all
set more off
set scheme sj
capture log close

// define user
global USER "`c(username)'"
display "Analysis run by $USER at `c(current_date)', `c(current_time)'"

// define working directory
global WDIR "YOUR WORKING DIRECTORY"

// create sub folders
capture cd "$WDIR/data"
if _rc! = 0 {
	mkdir "$WDIR/data"
}
global DATADIR "$WDIR/data"
store Latinobarometro 2017 data and catholics_LAC.csv into data sub folder

capture cd "$WDIR/tables and figures"
if _rc! = 0 {
	mkdir "$WDIR/tables and figures"
}
global OUTDIR "$WDIR/tables and figures"

// install required user packages
capture which tabout
if _rc {
	ssc install estout
}

// set working directory
cd "$WDIR"


**# 1. load data set
use "$DATADIR/Latinobarometro2017Eng_v20180117.dta"


**# 2. data preparation for independent variables
// mean center age
sum edad, meanonly
gen age_c = edad - r(mean)
label var age_c "Age"

// confidence in church
recode P14ST_C (1=3) (2=2) (3=1) (4=0) (-2=.a) (-1=.b), gen(conf_chrch)
label define church_lb 0 "no confidence at all" 1 "little confidence" 2 "some confidence" 3 "a lot of confidence" .a "NA" .b "DK"
label val conf_chrch church_lb
tab P14ST_C conf_chrch, m
label var conf_chrch "Trust in church"

// rating of pope francis
gen pope_francis = P61ST_3
mvdecode pope_francis, mv(-6 = .a \ -2 = .b \ -1 = .b)
label define pope_francis_lb 0 "very bad" 10 "very good" .a "do not know enough" .b "NA/DK"
label val pope_francis pope_francis_lb
tab P61ST_3 pope_francis, m
label var pope_francis "Rating Pope Francis"

// confidence in parliament, government and parties
foreach var of varlist P14ST_D P14ST_E P14ST_G {
	recode `var' (1=3) (2=2) (3=1) (4=0) (-2=.a) (-1=.b), gen(`var'_rec)
	label val `var'_rec church_lb
}
rename(P14ST_D_rec P14ST_E_rec P14ST_G_rec) (conf_parl conf_gov conf_parties)

// religious committment
recode S9_A (4=0) (3=1) (2=2) (1=3) (-3=0) (-2=.a) (-1=.b), gen(relig_commit)
label define relig_lb 0 "not practicing at all" 1 "not very practising" 2 "practising" 3 "very practising" .a "NA" .b "DK"
label val relig_commit relig_lb
tab S9_A relig_commit, m 
label var relig_commit "Religiosity"

// religious denomination
recode S9 (1=2) (2/5=3) (6/96=4) (97=1) (-2=.a) (-1=.b), gen(religion)
label define religion_lb 1 "None" 2 "Catholic" 3 "Evangelic" 4 "Other" .a "NA" .b "DK"
label values religion religion_lb
tab religion, m

// race/ethnicity
recode S10 (1=5) (2=1) (3=2) (4=3) (5=5) (6=4) (7=5) (-2=.a) (-1=.b), gen(race)
lab define race_lb 1 "Black" 2 "Indigenous" 3 "Mestizo" 4 "White" 5 "Other" .a "NA" .b "DK"
lab val race race_lb
tab race, m
label var race "Ethnic group"

// gender
gen female = sexo - 1
label define female_lb 0 "male" 1 "female"
label val female female_lb
label var female "Gender"

// education
rename REEDUC_1 education
label var education "Education"

// marital status (1=married; 0=single/separated/divorced/widowed)
recode S6 (1 = 1) (2/3 = 0) (-1 = .a), gen(married)
label define married_lb 1 "married/partner" 0 "single/separated" .a "DK"
label val married married_lb
label var married "Marital status"

// left-right self position
recode P19STC (-2 = .a) (-1 = .b) (-8 = .c), gen(leftright)
label define lr_lb .a "NA" .b "DK" .c "none"
label val leftright lr_lb
label var leftright "Left-right position"

// check all IVs
tab1 religion race female education married, m


**# 3. data preparation for dependent variable
// anthropogenic climate change denial
mvdecode P53N_I, mv(-1 = .a)
sum P53N_I
gen anth_climate_denial = r(max) + 1 - P53N_I
label define climate_anth_lb 4 "strongly agree" 3 "agree" 2 "disagree" 1 "very disagree" .a "DK/NA"
label val anth_climate_denial climate_anth_lb
label var anth_climate_denial "Anthropogenic climate change denial"
tab P53N_I anth_climate_denial, m


**# 4. empirical analysis
// generate country IDs
egen country_id = group(idenpa)

// define sample
qui reg anth_climate_denial i.conf_chrch pope_francis i.female c.age_c ib(5).education i.race i.married c.leftright i.relig_commit ib(3).religion
keep if e(sample)

// descriptive statistics 
gen year = 2017
gen date = mdy(mesreal, diareal, year)
format date %td
bysort country: gen respondents = _N

**# 4.1 Table 1
estpost tabstat edad leftright pope_francis, stats(min median mean max sd) columns(statistics)
esttab using "$OUTDIR/table_1.rtf", cells("min(fmt(2)) p50 mean max sd") nomtitle nonumber replace

foreach var of varlist anth_climate_denial conf_chrch education relig_commit race married religion {
	estpost tabulate `var', nototal
	esttab using "$OUTDIR/table_1_`var'.rtf", cells("pct(fmt(2))") nonumber nomtitle noobs replace
}

// (multivariate) models
global CONTROLS i.female c.age_c ib(5).education i.race i.married c.leftright i.relig_commit

qui {
	meologit anth_climate_denial ib(2).religion i.conf_chrch pope_francis || country_id:
	est store model_1

	meologit anth_climate_denial $CONTROLS || country_id:
	est store model_2

	meologit anth_climate_denial ib(2).religion i.conf_chrch pope_francis $CONTROLS || country_id:
	est store model_3
}

// omnibus F-test trust in chruch
est restore model_3
test (1.conf_chrch=0) (2.conf_chrch=0) (3.conf_chrch=0)

// interactive model
meologit anth_climate_denial ib(2).religion c.conf_chrch##c.pope_francis $CONTROLS || country_id:
est store model_4
	
// testing for significantly different effects in Argentina
gen arg = country_id == 1

meologit anth_climate_denial ib(2).religion c.conf_chrch##c.pope_francis##i.arg $CONTROLS || country_id:
est store model_5


**# 5. export tables and figures
**# 5.1 Supplementary Table 1
preserve
collapse (min) date_min = date (max) date_max = date (median) respondents, by(country_id)
list
restore

**# 5.2 Supplementary Table 2
esttab model_1 model_2 model_3 model_4 model_5 using "$OUTDIR/Supplementary_Table_2.rtf", replace label ///
						cells("b(fmt(3) star) se") star(+ 0.1 * 0.05 ** 0.01 *** 0.001)  ///
						mtitles("Model 1" "Model 2" "Model 3" "Model 4" "Model 5") ///
						order(*religion *conf_chrch pope_francis *relig_commit leftright age female married race education)

**# 5.3 data for figure 2				
// retrieve and export model coefficients
est restore model_3
mat b = e(b)
mat se = vecdiag(e(V))
mat export = b', se'

preserve
clear
svmat export
local names: rownames export
gen var = ""
forval i = 1/`: word count `names'' {
   replace var = `"`: word `i' of `names''"' in `i'
}
rename export1 beta
rename export2 se
replace se = se ^.5
export delimited "$OUTDIR/model_coef.csv", replace
restore

**# 5.4 data for figure 3
// retrieve and export predicted probabilities
est restore model_3
foreach x in 2 4 {
	margins, predict(fixedonly outcome(`x')) at(female = 1 education = 5 relig_commit = 2 race = 3 married = 1 leftright = 5 pope_francis = 5 religion = (2 3))
	mat res_rel`x' = r(b)', vecdiag(r(V))'

	margins conf_chrch, predict(fixedonly outcome(`x')) at(female = 1 education = 5 relig_commit = 2 race = 3 married = 1 leftright = 5 pope_francis = 5)
	mat res_`x' = r(b)', vecdiag(r(V))'
	
	margins, predict(fixedonly outcome(`x')) at(female = 1 education = 5 relig_commit = 2 race = 3 married = 1 leftright = 5 pope_francis = (0 10))
	mat res_pope_`x' = r(b)', vecdiag(r(V))'
}
mat res = res_rel2 \ res_rel4 \ res_2 \ res_4 \ res_pope_2 \ res_pope_4
preserve
clear
svmat res
rename (res1 res2) (beta se)
replace se = se^.5
export delimited "$OUTDIR/model_pred_prob.csv", replace
restore

**# 5.5 data for figure 4
// retrieve and export conditional predicted probabilities
meologit anth_climate_denial ib(2).religion i.conf_chrch##c.pope_francis $CONTROLS || country_id:

margins, predict(fixedonly outcome(4)) at(conf_chrch = 3 female = 1 education = 5 relig_commit = 2 race = 3 married = 1 leftright = 5 pope_francis = (10(-1)0))
mat res_int = r(b)', vecdiag(r(V))'

preserve
clear
svmat res_int
rename (res_int1 res_int2) (beta se)
replace se = se^.5
export delimited "$OUTDIR/model_pred_prob_int.csv", replace
restore

								
**# 6. jackknife leave-one-out
forval x = 1/18 {
	qui meologit anth_climate_denial ib(2).religion c.conf_chrch##c.pope_francis $CONTROLS if country_id != `x' || country_id:
	est store model_6_`x'
}

**# 6.1 Supplementary Table 3	
esttab model_6_* using "$OUTDIR/Supplementary Table 3.rtf", replace label ///
						cells("b(fmt(3) star) se") star(+ 0.1 * 0.05 ** 0.01 *** 0.001)  ///
						order(*religion *conf_chrch pope_francis) ///
						keep(*religion *conf_chrch pope_francis c.conf_chrch*)